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Abstract — The magnetotelluric (MT) method has become more 
widely used in hydrocarbon exploration. The inversion of MT 
data, which can determine the electrical structure of subsurface, 
is a nonlinear and multimodal optimization problem. Particle 
swarm optimization (PSO) algorithm is a good solver for this 
geophysical inversion problem, whereas it has a shortage of 
heavy computation time. An improved parallel adaptive PSO 
inversion algorithm for MT data is proposed in order to decrease 
the computation time. The performance of the proposed 
algorithm was evaluated on the Dawn 4000L supercomputer 
using the synthetic MT data of ID layered geo-electrical models 
of three and four layers. The numeric results show that the 
proposed algorithm can obtain the as good solution as the serious 
PSO inversion algorithm, and can reduce the computation time 
obviously when more computing nodes been employed. This 
result indicates that proposed improved parallel inversion 
algorithm can deal with the computation time problem and 
provide theory and technology support the MT data non-linear 
inversion based on PSO. 
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I. Introduction 

The magnetotelluric (MT) method utilizes naturally 
occurring low-frequency electro-magnetic waves to determine 
the electrical resistivity of the Earth's subsurface. A knowledge 
of the resistivity is important since rock types important to 
hydrocarbon exploration can be differentiated on the basis of 
resistivity value [1, 2]. It has some advantages such as low cost, 
large probing depth, not affected by high-resistivity shielding 
and high resolution to the low-resistivity layer. In the last 
decades, the MT method has become more widely used in 
hydrocarbon exploration [3]. 

Data inversion is one of the core issues of the MT method, 
which is classified into linear and nonlinear categories. To 
overcome the flaws of the classical linear inversion method, 
some nonlinear inversion methods, such as Monte Carlo, 
simulated annealing and genetic algorithm are employed to the 
MT inversion [4]. However, they have some disadvantages [2]. 

The particle swarm optimization (PSO) is a relative novel 
global optimization algorithm based on swarm intelligence [5]. 
It has several advantages including fast convergence, strong 
capability of global optimization and simple parameter 


adjustment, and has been successfully applied to many fields 
[6]. In recent years it has also been introduced into the 
geophysical inversion [7-11]. Although the disadvantage, such 
as local optimum, of geophysical inversion based on PSO has 
been improved by employing the adaptive inertial weight, it 
still has the problem of heavy computation time when the 
forward modeling becomes more complex [2]. Parallel 
implementation can reduce the computation time of inversion 
obviously. 

Message Passing Interface (MPI) is one of the most popular 
parallel computing software environment [12]. An 
asynchronous parallel adaptive PSO inversion algorithm based 
on MPI for the MT oil-gas exploration data is proposed in this 
paper. It is applied to the inversion of the synthetic MT data 
and its parallel performance is analyzed. 

II. Methodology 
A. Forward modeling 

Assume an ID layered model through a cross section of n 
layers, of which from up to down the resistivities are 

p v p 2 ,..., p n , respectively, and the depths are 1\ , , . . . , h n , 
respectively, where h n = oo . For such a ID layered model, 

apparent resistivity p a and phase (j) a can be calculated as 
follows: 


PjO)) = 


\ Z(co) I 2 


(j) a = arctan 
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Where CO — 2 nT is angular frequency, p is magnetic 
permeability, Z(co ) is wave impedance on the surface that can 
be calculated by the following recurrence formula: 
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Where k j — ^ Jicop / p { is the complex wave numbers of 
the z-th layer, Z Qi is the characteristic impedance of the z-th 
layer, and Z. is the wave impedance of the top of the i-th layer. 
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The presentation above indicates that apparent resistivity 
and phase on N-layer geo-electrical section can be expressed 
by functions of signal periods and section parameters: 

Pa ~ ( 3 ) 

(fra — fl ’ Pi ’•••’ Pn ’ ’ ^2 ’ • • •’ h N - 1 ? T} 

B. Improved PSO 

In the PSO algorithm [5], the optimal solution is obtained 
through collaboration among individual which represented by 
particles. Suppose that the search space is of n dimensions, 
total number of particles is m, the position of the i - th particle in 
the n-dimensional space is xi, and its flight speed is vi. Each 
particle has an adaptive value determined by an optimization 
target function, and is aware of the best position pi that is found 
till now by itself and its current position, as well as the optimal 
position pg of the whole swarm found till present. Then the 
speed and the position of each particle can be updated by the 
following formulas: 

vf +1 =wvf +c 1 r 1 (p i —x*) + c 2 r 2 (p g -xf) 

(5) 

•y.fc+1 + l (fsA 

X i ~ X i +V / I 6 ! 

k 

where V- is the flight velocity of the i - th particle after k-th 

k 

iteration, x. is the position of the i - th particle after k-th 
iteration, r x and r 2 are random parameters evenly distributed 

in 0~1, c t and c 2 are weight factors, and wis inertial weight 

which is important to the balance between global exploration 
and local exploitation. 

The Adaptive PSO [2] adjust the inertial weight 
dynamically according to the particle velocity of the population 
as follows: 

if v L g > v e > then w (k + 1) = w(k) / p ; 

if V Lg <V e ’ then w(k+ \) = w(k) * p ; 

(7) 

if V avg = V e ’ then W (k + 1) = w(k) , 

where 

m n 

v l g =SZ lv S | ) / ( m *«) 

i = 1 j = 1 

is the average velocity of the k-th generation population, 

v ^-(2 k/iT^-TO) 2 

is the expected velocity of the &-th generation population 
which attenuates exponentially, here T max is the maximum 

evolution iteration, 7j is the local exploitation iteration, 


T max — T x the iteration of global exploration, and parameter p 
is the change rate of inertial weight. 

C. Parallel Strategy 

Some results [12] reported show that the coarse granularity 
is more effective then the fine granularity, because it can 
reduce the communication time so as to gain better 
performance. So we choose the coarse granularity in our 
parallel implementation of inversion, that is to say we assign a 
number of particles to one processor. 

We regard the particles in the same processor as an 
independent sub- swarm, which has the full information. The 
particles in the same processor (sub-swarm) can construct the 
optimal solution independently guided by their own 
information. Different sub- swarm (different nodes) interact 
each other by exchanging the information such as sub-swarm 
optimal solution. We name the dividing strategy as coarse 
granularity interacting multi particle swarms. 

In order to reduce the communication time of parallel 
computing, we choose the asynchronous parallel strategy, 
which means each node (sub-swarm) run certain iterations 
locally and exchange the information each other controlled by 
the master node and then begin the next certain iterations. 

D. Improved Parallel PSO inversion algorithm 

The forward model can be written as d = A(m) , where 

m is a model parameter, A is forward functional, and d is the 
theoretical value corresponding to the model m. Inversion is to 

solve the model parameter m from the observed value d obs , 
which makes the fitness error between the theoretical 

value d = A(rri) and the observed value d obs least. The 
objective function of inversion is defined as the norm L 2 of the 
difference between the observed and theoretical values to 
describe the fitness degree, which is expressed as 

P(m ) =11 d obs - Aim) II 2 min 

( 8 ) 

For an N-layer geo-electrical section, the model parameter 
m is a (2N-1) dimensional vector 
(p v p 2 ,...,p N ,h [ ,h 2 ,...,h N _ l ) T . If the observation is 

conducted at K frequencies, the observed value d° bs will is a 
2* K dimensional vector (p x , p 2 ,...p K ,(f> l ,(f} 2 ,...(f} k ) r , where 

p { , (fr i are apparent resistivity and impedance phase at the i-th 

frequency. The inversion of MT data for ID N-layer model is 
to find the model parameter 

m = (p x , p 2 ,..., p N ,h x ,h 2 ,...,h N _ l ) T , which results a best 

match between the calculated value d = A(jri) from Eqs. (3) 

and (4) and the observed value d ° bs . 

The improved parallel PSO inversion algorithm is described 
as follow: 
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Step 1: Input node number (sub-swarm number) N, particle 
number S of each sub-swarm; 

Step 2: Master node generate N*S parameters of initial model 
randomly, corresponding to initial positions of N*S 
particles, and send z'-th sub-swarm to z-th node; 

Step 3: Each node calculates each particle’s theoretical value 
Aim ) of their own sub-swarm respectively, using 
formulas (3) and (4); 

Step 4: Each node calculates objective function value of each 
particle of their own sub-swarm by formula ( 8 ); 

Step 5: Each node selects the best particle of their own sub- 
swarm; 

Step 6: Each node calculate current average velocity V avg of 

sub-swarm, and adjust inertial weight w of sub-swarm 
dynamically according to v avg and v e by formula (7); 
Step 7: Each node adjusts the position and velocity of each 
particle of their own sub-swarm, according to 
following formulas (5) and ( 6 ); 

Step 8: If every nodes runs certain iterations, they exchange 
the information such as the position and velocity of the 
best particle of each sub-swarm, controlled by the 
master node; 

Step 9: if global convergence or maximum number of 
iteration is met, go to Step 10, otherwise return to Step 
3 to perform next iteration; 

Step 10: Master node output result of inversion, finish 
calculation. 

III. Numerical result and analysis 

The parallel PSO inversion algorithm is conducted to the 
noise-free MT data, with 20 particles, maximum iteration 
number 2000 and learning factor cl=c2=2. In inversion, the 

value ranges taken are p x - 0~500Qm, P 2 - 0~100Qm, p 3 = 

0~4000Qm, /?i=0~ 2000m, /z 2 =0~4000m, respectively. The 
inversion result is shown in Table 1, where the sum of relative 
errors is that of inversion model parameters 

(p x , p 2 , p 3 ,h x ,h 2 ) and real model parameters. 

MT data inversion by parallel PSO is made on the four- 
layer (type HH) model. The result is listed in Table 2. In 
inversion, 20 particles, maximum iteration number 2000 and 
learning factor cl=c 2=2 are adopted, and the value ranges 

taken are =0~1000Qm, /) 2 = 0 ~ 1 OOOQm, p 3 = 0~1000Qm, 

P\ - 0~1 OOOQm, /?i=0~4000m, /z 2 =0~4000m, /? 3 =0~ 4000m, 
respectively. In these ranges, 20 initial particles are generated 
randomly for iterative inversion. 

Tables 1 and 2 show the accuracy of result of parallel 
inversion is as well as that of serial one’s but it spends much 
less computation time. 

We plot the speedup and speedup efficiency of parallel 
PSO inversion algorithm in figure 2 following the data of 
Table 2. From Fig. 2, we can find that more computation 


nodes lead to less computation time but less speedup 
efficiency, for inversion. 


TABLE I. Parallel inversion results on the three-layer model 


Model 

Parameter 

Real 

Mode 

1 

MT Inversion 

1 

node 

10 

nodes 

20 

nodes 

30 

nodes 

pi 

110 

110.1 

110.0 

109.9 

109.7 

p 2 

20 

20.01 

20.0 

20.43 

20.1 

p3 

1200 

1199.8 

1200.0 

1200.2 

1199.8 

hi 

500 

498.6 

498.9 

498.2 

500.2 

h 2 

2000 

2001.3 

2001.6 

2015.7 

2017.3 

Running 

number 

- 

50 

50 

50 

50 

Calculation 

time(s) 

- 

482.9 

129.4 

90.9 

75.3 


TABLE II. Parallel inversion results on the FO UR-layer model 


Model 

Parameter 

Real 

Mode 

1 

MT Inversion 

1 

node 

10 

nodes 

20 

nodes 

30 

nodes 

pi 

100 

99.9 

100.2 

99.9 

99.9 

p 2 

20 

19.2 

20.6 

20.0 

19.9 

p3 

300 

281.6 

300.4 

286.5 

319.4 

P4 

10 

10.0 

10.4 

9.9 

10.1 

hi 

600 

612.1 

595.9 

599.5 

601.4 

h 2 

1500 

1468.4 

1502.6 

1537.2 

1512.1 

h3 

3000 

3176.9 

2899.2 

3018.3 

2979.4 

Running 

number 

- 

50 

50 

50 

50 

Calculation 

time(s) 

- 

617 

132.2 

93.8 

79.8 


IV. Conclusion 

The MT method has become more widely used in 
hydrocarbon exploration. In order to avoid the disadvantage of 
heavy computation cost, a synchronous parallel adaptive PSO 
inversion algorithm is proposed based on MPI. The 
performance of the proposed algorithm was evaluated on the 
Dawn 4000L supercomputer using the synthetic MT data of 
ID layered geo-electrical models of three and four layers. 

The numeric experiments shows that the proposed parallel 
inversion algorithm can obtain the as good inversion solution 
as serious version, and can reduce the computation time 
obviously when more computing nodes been used. This result 
indicates that proposed asynchronous parallel inversion 
algorithm can deal with the computation time problem and 
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provide theory and technology support the MT data non-linear 
inversion based on PSO. 
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Figure 1 . THE SPEEDUP AND SPEEDUP EFFICIENCY OF FOUR 
LAYER MT DATA PARALLEL INVERSION. 
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